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ABSTRACT 


Six vertical discretization schemes are compared for a 
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Variables are defined for two staggered and one unstaggered 
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results of the finite element models are not more accurate 
than the corresponding finite difference results.  Oscilla- 
tions occur in the temperature profiles near the surface for 
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I. INTRODUCTION 


The finite elements method (FEM) has been used in several 
atmospheric prediction models (Staniforth and Mitchell, 1977; 
Williams and Schoenstadt, 1980, Beland, et al., 1983). The 
FEM is a Special case of the Galerkin procedure in which the 
dependent variables are approximated by a finite sum of 
spatially varying basis functions with time dependent coeffi- 
cients. The FEM basis functions are low order polynomials 
which are zero except in a localized region. The Galerkin 
procedure produces a set of coupled ordinary differential 
equations which are solved by using finite differences in 
time (Haltiner and Williams, 1980). 

FEM models are potentially more accurate than finite 
difference method (FDM) models. Winninghoff (1968), Arakawa 
and Lamb (1977), and Schoenstadt (1980) have demonstrated 
the advantages of spatial staggering of predictive variables 
in finite difference approximations of shallow water equa- 
tions. The results of Williams and Schoenstadt (1980) indi- 
cated that FEM models should either use staggered nodal points 
in the momentum equations or’ unstaggered nodal points in the 
vorticity-divergence equations. 

The purpose of this research is to compare six linear, 
baroclinic, vorticity-divergence equation models. A finite 


difference and a finite element model are written for each of 


1 2 


the three vertical grids depicted in Figure l. Grid A, which 
was originally developed by Lorenz (1960) for energy conser- 
vation, 1S a widely used grid for finite difference models. 
However, Tokioka (1978) has shown that this grid has a 
computational mode in the temperature field, and Arakawa (1984) 
has found a false small scale baroclinic instability with this 
grid. Tokioka (1978) also found that grid C has computational 
modes in all fields. Béland, et al., (1983), analyzed the 
finite element formulation for grid C and found noise gener- 
ated by certain forms of friction. Grid B was introduced by 
@marney and Phillips (1953) for a finite difference quasi- 
geostrophic model.  Tokioka's (1978) analysis has shown that 
this grid has no computational modes. It is hoped that the 
finite element formulation for this grid will have high 
accuracy without the problems the other grids have had with 
computational modes. 

Each of the models used in this thesis are derived from 
the governing equations described in Chapter II. The models 
are one-wave spectral in the horizontal and have a fixed lid 
at the top of the atmosphere. Diabatic heating and mountain 
topography can be included in each model experiment. 

The results of three experiments with analytic initial 
conditions will be examined. As described in Chapter III, the 
experiments are 1) an initial perturbation in the meridional 
flow, 2) flow over mountains, and 3) flow with a mid-atmosphere 


diabatic heat source. Each experiment will be presented with 
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GRID C 


Three vertical grids to be compared for N- 
layer models. Perturbation variables are 
defined at either unstaggered levels (solid 
lines at height Z2') or staggered levels 
(dashed lines at height Z). C is vorticity, 
D is divergence, T is potential temperature, 

T is basic state potential temperature, 

$ is geopotential, 9s is surface geopotential, 
u is east-west velocity, u is basic state 
east-west velocity, v is north-south velocity, 
Z is vertical velocity, Q is dibatic heating, 
and MTS is forced vertical velocity. 
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six and sixty vertical layers. The results of the 60-layer 
models will show if the models are converging to the same 
solution. For each model, comparison of the six and 60- 
layer results indicate how well the lower resolution model 
approximates the convergent solution. The purpose of these 
experiments is to identify any problem areas in the models 


and to demonstrate the general characteristics of each model. 
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ll. MODEL DESCREBPEBISNS 


A. MODEL FEATURES 

The six numerical models are developed with several fea- 
tures which make the models easy to modify for a wide range 
of experiments. The depth of the atmosphere, number of 
levels and the depth of each layer are variable in each model. 
Diabatic heating and forced vertical velocity due to mountain 
topography are included in the governing equations. The user 
can prescribe heating or forced vertical velocity functions; 
or run the model with these terms defined as zero. The models 
are written in modular structure using FORTRAN '77. There is 
parallel construction between models. The subroutines used 
in one model are very similar to those used in the other five 
models. The models run very quickly on an IBM-3033 mainframe. 
A 96-hour forecast for a 12-layer FEM uses less than five 


seconds of computer processing time. 


B. GOVERNING EQUATIONS 

Each model approximates the same set of governing equa- 
tions. The vorticity (2-1), divergence (2-2), surface geo- 
potential (2-3) equations and the first law of thermodynamics 
(2-4) are the prognostic equations for the forecast variables 
vorticity, divergence, surface geopotential and potential 
temperature. The surface geopotential equation is the lower 


boundary condition on the vertical velocity. The vertical 
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coordinate Z = -ln(p/Po) is used, but the non-Boussinesq 
terms involving ae are replaced by one. The prognostic 


equations in the coordinates x,y,Z,t are 


dé ομως S 
me DF BV toe ὃς ὃν 3y "et ο 
dD ο 2 αν 2 du ὂν ΟΖ du 
το * ‘ay’ t? 3. 3x * 3x 3z 
97 ov 2 " P 
ο στ τσ D a 
ας. 2-3 
Em - BI, (2 3) 
S = 
D °. E 
Here, 
M. ον ou 
ζ is the vertical component of vorticity, ¢ = Πα ὃν’ 
D is the horizontal divergence, D = E + A 


@ is the geopotential, ġġ = gZ, 

js 15 the surface geopotential, 

T is the potential temperature, 

u is the x-component of velocity, 

V is the y-component of velocity, 

Z is the vertical velocity, 

Q is the diabatic heating per unit time per unit mass, 


MTS is the forced vertical velocity due to flow over 
mountain topography, discussed in Chapter III, 
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is dt 
dy” 


|> vU Fh 


Qu 


t - 
v2 


9 
dy 


9 9 ; 
τ πα. + 2 


is the Coriolis parameter, 


= 
3n 


is the horizontal Laplacian operator. 


The prognostic equations are linearized by expanding the 


variables into the following mean state and perturbations: 


lee te) = 


Vix. 2, 6) 


Z(x,Z,t) 

πο ο ε) 
o(x,y,Z,t) 
o,(x,y,t) 


Cx t) M 


Dix Zt) 


O (2, 2 E) 


MTS (X€) 


u(Z) + u'(x,Z,t) 


Vox) 


Z(x,Z,t) , 


τίν,α) : ΠΠ 


$032) +o (k Ao 


, 


Q ° (y) + jo (x E) 


C (X AmE) 


tx. ut) 


O x 2 t) 


MTS' (xD) 


1. 


, 


(2-5) 


(2:299) 


(2:599 


(2-8) 


(2290) 


(25 


(251. 


(2-1 


(2-15 


(2-14) 





The vorticity and divergence in this system are ζ 3x' 
t 
D = Lm The linearized forecast equations are 
oc! = iui a. = ΕΝ ES t = 
TEN fD Ux p vs οσο) 
PE ρε - u 2D! sy gun - 22 We οφ: (2-16) 
ο Š ox 3X dZ πας, 
ðX 
um — o V NUT 2 OT Ν 
TM EN ab, .. 
DES = = > ULT 3x SUM Dm =- RTZ + MTS , (2-13) 


where R is the gas constant for air. 

The diagnostic variables, u', v', 2', $', are calculated 
from the forecast variables using the definitions of divergence, 
vorticity, the hydrostatic equation and the continuity equa- 


tion. The relations are 


s penc 220502 23:9) 

= TF (2-20) 

oo = RT, (2-21) 
x 

D' + at = 0. (2= 22) 


The use of primes to denote perturbation quantities will be 


discontinued. All quantities used in the remainder of 
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the paper will be perturbation quantities unless otherwise 
noted. 

The mean state is assumed to be in hydrostatic and geo- 
strophic balance. The term 93T/93y in the first law of thermo- 
dynamics can be evaluated by taking 3/3y of the hydrostatic 


equation and substituting for 3$/98y from the geostrophic 


relation, 9$6/3y - -fu. Thus, 
T f əu 
Sy | Re Ll 


Geostrophic balance of the mean state at the surface implies 


T» H 
ὃν = -= fu cc ; (2-24) 
The expressions (2-23) and (2-24) are substituted into equa- 


tions (2-17) and (2-18), respectively. 


A singlewave spectral representation is used in the 


x-direction, with wave number u = 21/L, where L is the wave- 
length in the x-direction. The perturbation quantities have 
the form 
BVA Kee δ- A, (Z,t) cos ux + A,(Z,t)sin ος ; (2-25) 
De oe). = Di(2,t)cos ux + 522, Έ) 51η Wore (2-26) 
E AE S T, (2,t)cos ux + T,(Z2,t)sin lc (2-279 
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$2 (x, t) = S, (t)cos ux + S, (t)sin με , (2-28) 


m(x,2,t) = U, (Z,t)cos ux + U,(Z,t)sin μετα, (2-29) 
EX 2,t) = V4 (2,t)cos ux + V5 (2,t)sin ucc (2-30) 
E o = W, (2,t) cos ux + W5(2,t)sin UX , (2-31) 
ο σσ, ο) = Hi (2,t)cos ux + H,(Z,t)sin "ασ. (2-32) 
mex, 4,t) = Q4 (2, t)cos ux + Q,(Z,t)sin Uc c2 99) 
MTS'(x,t) = MTS, (t)cos ux + MTS, (t)sin NX . (2-34) 


The relations (2-25)-(2-34) are substituted into equations 
(2-15)-(2-22). The prognostic and diagnostic equations are 
separated into equations for the cosine and sine terms. The 


resultant prognostic equations are 


B . "E _ - 
3A, ü 
aD — 

m _ τσ _ "au 2 °. 
E A) στ uuD, BU, ~ kag Wo + u H) , (2-37) 
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2E -- _ du 2 
oT m m 
να. f du _ 9T 
3p ο Tui t RaZ 10 32850 EN ον 
aT T = 
υπο f du _ oT 
jt 7 U"U"Lus nan 2 RM R E 
951 _ _ - 
στ = TuS, * fuV, - RTW, + MTS, , (2-409 
952 m αι - 
The resultant diagnostic equations for u and v are 
D 
ea eS 2 
Ui = mt f (2-43) 
D 
=e = 
U, πό (2-44) 
vo «ΞΡ se; (2-45) 
1 ΠΠ. 
A 
a fae a 
V, = k (2-46) 


Geopotential values above the surface are obtained by Sima 
grating the hydrostatic equation from the surface (Z = 2ο) to 


height 2: 


2 


H = R Jj Diz eazy (2-47) 
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Z 


D MM uos, (2-48) 
Z 
0 
The vertical velocity is calculated by integrating the con- 
tinuity equation from the top of the atmosphere (Z = Zn) down 
to height Z. The upper boundary condition, Z = 0 at Z = Zr 


is used, This boundary condition is not exact, but some form 
of it is used in most numerical models. The diagnostic 


equations for the vertical velocity are 


ap 
Ww, = J D, (Z,t)dZ , (2-49) 
an 
Wy = J D,(Z,t)dzZ . (2-50) 
Eguations (2-35)-(2-50) are the prognostic and diagnostic 
equations that govern all six numerical models. Using the given 


basic state and the one-wave spectral perturbation quantities, 
the governing equations reduce to functions of Z and t. The 
models are effectively one-dimensional. 

To display the results of each model, the sine and cosine 
amplitudes of each variable are combined to determine the 
amplitude and phase of a single cosine wave in the x-direction. 


A typical variable has the form 


Wee ©) RUASME)COSTQUNx-OÓ) , (2-51) 
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where the amplitude is A(Z,t) and the phase is 6(Z,t). The 
amplitude and phase are calculated at each level for all 


variables. 


C. TIME DIFFERENCING 

Two forward time steps are taken to start each model and 
then leapfrog time differencing is used. The leapfrog scheme 
is employed because of its ease to code. A Robert filter is 
used to reduce the amplitude of the computational mode gener- 
ated by the leapfrog time differencing. The filter is dis- 
cussed by Haltiner and Williams, 1980. For a prognostic 
variable F, calculate Βλ, the average value of F at time 
step (n-l)At, using equation (2-52). Using the unaveraged 
values at time step nât, compute the tendency (s= . from aes 
predictive equation. The predicted value at time step 


(n+l) At is calculated using equation (2-53). 


Eoo = Fs) POE HS NM (2-52) 
where Υ 15 ἃ weighting function. 
F = F _ + 2At(2Ë (2-53) 
τ πι ος ῃ 
In all thesis experiments, y = 0.05 is used. The time step 


for each experiment is calculated in the model by requiring, 


for computational stability, 
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våt = š , (2-54) 
where v = uc, and c is the typical phase speed of an external 
gravity wave. 

D. VERTICAL GRIDS 
Each of the models uses one of three vertical grids. The 


v 


three ways of distributing the variables over discrete levels 
are depicted in Figure l1. The staggered levels are represented 
by the dashed lines in Figure 1. Notice that the heights at 
which the variables are defined change between the three 
grids. The notation used in this paper to denote the staggered 
and unstaggered levels is consistent with the conventions used 
in the coded models. The height of the unstaggered levels 
is denoted as Z'. The height of the staggered levels is 
denoted as Z. In the models, both Z4 and Σι are defined to 
be the surface of the earth. It is assumed that the staggered 
Never αι is exactly in the middle of the layer between HY 
and Zi. This distinction is important because the models can 
have layers with unequal depth. Thus, the height of the 
Staggered levels is defined relative to the height of the 
unstaggered levels. 

A finite difference model is written for each of the grid 
structures. The models are denoted as DFM-A, FDM-B and 


FDM-C. Similarly, finite element models using the three grids 


are indicated by FEM-A, FEM-B and FEM-C. 
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E. FINITE DIFFERENCE MODELS 

The only differences in the equations between the three 
FDM models are the approximations of terms involving du/dzZ 
and IT/9Z in the prognostic equations and the approximations 
of the integral in the diagnostic geopotential equation. 
Centered difference approximations are used, except at the 
boundaries where one-sided differences are employed. The 
finite difference approximations used in the prognostic equa- 


tions are listed in Appendix A. 


F. FINITE EDGCEMENT MOBEDS 
l.  FEM-EÉ 
The unstaggered FEM model is the simplest of the three 
FEM models. Each of the dependent variables is expanded into 
a finite series in terms of the eigenfunctions $5 (2). The 
eigenfunctions for this model are depicted in Figure 2. The 


eigenfunction expansion for a typical term is 


NFI 


š j a 
A,(Z,t) = l Ai (€) $4 (2) : (2255p 


J=1 


The finite element approximations for the vorticity, 
divergence and thermodynamic equations are derived by substi- 
tuting the eigenfunction expansion for each dependent variable 
into equations (2-35)-(2-40). Each equation is multiplied 


by $1 (2) and integrated with respect to Z from the bottom to 
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Figure 2. Eigenfunctions defined at unstaggered levels 


(solid lines) for grid C. Dashed lines 
represent staggered levels. 


the top of the atmosphere. Each term in the equations is the 
finite sum of separate integrals. Only the integrals of 
overlapping eigenfunctions are nonzero. The resultant equa- 
tions, listed in Appendix B, are matrix equations. For an 
N-layer model, the variables are calculated at N+l discrete 
levels and the variable can be represented as a column vector 
of length N+l. Thus, the forecast equations for vorticity, 
divergence and potential temperature become vector equations 
which contain four (N+l) x (N+1) matrices multiplying appro- 


priate vectors. 


2 7 


The mass matrix for this model is defined by 


ΠΣ J $,020,(2)d2, for i - 1,...,N4l 


The matrix N is defined for terms multiplied by u, 


v 


En EN ë 
N = u. 6. (Z) 6, (2) >. (2) dz, 
τ j=i-1 k=í-1 Jz ο σα, 


for i= e. nil! 


The matrix P is defined for terms multiplied by du/dZ, 


a ad. 
r s | u, [f apt Oy (2) 9, (242, 


for i z olye. NPI 


The matrix R is defined for terms multiplied by 91/97, 


ER Y _ E ° 
R = | T: Ó (σι ο (2) g2 
τ πο πρ -. πας, 


for 1i τ m 


The method to evaluate the terms of matrices (2-56)-(2-59) 


(2-567 


(22578 


(2=5 p 


(2-597 


and the resultant formulas for the matrix elements are listed 


in Appendix C. 
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The vorticity, divergence and thermodynamic equations, 
written in matrix (indicated by a double underline) and 


vector (indicated by a single underline) form are 


M-S=(AL) L G ο ο ο Ν.Α, (2-60) 
Μ. ασ) - M:(-£-D2-8-V2) + u-N-Al , (2-61) 
M-Se(DL) = M-(f-AL-8Ul) - wN-D2 - u:P-W2 
2 
+ u M:H1 , (2-62) 
M-S-(D2) = M+-(f:A2 -8U2) + u-N-Dl + u-P-W2 
2 a 
ο να, (2-63) 
d = f 
M-S-(Tl) = -u-N:T2 + (É) P-Vl - R-Wl + M-Ql , (2-64) 
M-S-(T2) = i-N-Tl & (D-P-V2 - R-W2 4 M-Q2 . (2-65) 


Equations (2-60)-(2-65) are simplified by multiplying 


each equation by M1 and applying the Robert filter. The 


i z1 = z1 
Maerices M +N, ο ee eee :εἱἵδπεᾶάτε5,  Ίπεν are con- 


Im 


structed in the initialization subroutine and stored for use 
in the forecast subroutine. The matrices are multiplied by 


the appropriate vectors with values for time level nAt. The 
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resultant forecast equations are vector equations and the fore- 
cast value for the i-th vertical level is the sum of values 

in the i-th location of each vector in the equation. The 
prognostic equations for the vo i divergence and poten- 


tlal temperature vectors are 





























rz as a A 
Al.,j' 7^ ἈΙΒΆΒ 1 + 24t(-f DIC orv m Meen NE (2-66) 
_ -1 
Αλ = A@BAR, , * 2at(-f:D2-&-V2 41-M "-N:AI) , (2767) 
= | f= ee 
Dl |, 7 DIBAR, , * 2At(f-Al- 8ü-Ul -u*M "N-D2 
- wx l-pw2 6 ug) , (2-68) 
D 
— ° ° “LR e 
D2 ., = D2BAR _, + 2At(£-A2 = 6°02 Fil Mion om 
-1 2 
πμ πη a , (2-69) 
Tl νη... ποι; 
==n+1 — =i: H m r ο 
€ poule vico ailes» o Qu, (2-70) 
. ον a 
T2 į] 7 T2BAR ο ο ο N 
E (=) -M "+P +2 = M ~ απο το... (2-71) 
m EE πο 














where the subscripts n+l, n, n-l, refer to the values of the 


vectors at time step (nt+l)At, nAt, and (n-1)At, respectively. 
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The surface geopotential and the diagnostic variables 
are calculated using the same equations that are used in 
the unstaggered model FDM-C. 
2. FEM-A 
The FEM-A model defines vertical velocity at the un- 
staggered levels in terms of the eigenfunctions 34 (2). The 
other variables are defined at the staggered levels in terms 


of the eigenfunctions $4 (2). The eigenfunctions for this 


model are depicted in Figure 3. 
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Flgure 3.  Eigenfunctions for grids A and B. Egienfunc- 
tions ψί(2) are defined for the unstaggered 
levels (solid lines at height Z') and 9$(2) 
are defined for the staggered levels 
(dashed lines at height Z). 
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The finite element approximations for the vorticity, 
divergence and thermodynamic equations are derived using the 
method described for the FEM-C model. The resultant Galerkin 
formulation of the equations, listed in Appendix D, are matrix 
equations. The mass matrix, M, and the matrix, N, for terms 
of the form u multiplied by variables defined at the nodal 
points óf $4 (2) have the same formulas as M and N in the FEM-C 
model, equations (2-56) and (2-57). However, the eigenfunc- 


tions $4 (2) are not defined at the same levels in both models. 


The matrix Q is defined for terms ot Che τ η soy, 
š ! Z 
v EN E f T d$. 
O Ξ i u. — ġ ġ:dZ . (2-72) 
- "ΡΨ, 
The matrix P is defined for εις ου σας amm oe, 
; Z 
Jill sr s A | T ae. 
P = ) yo rece 9 p,o.dZ . (2-73) 
= π᾿ 
The matrix R is defined for terms ΟΕ ΕΞ, ΗΠ SW, 
; l Z 
PEL το - p ai 
E a / 37 UQ0,d2 . (2-74) 


The staggered eigenfunctions present three general 


problems for evaluating the elements of the five matrices. 
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First, for an N-layer model, portions of eigenfunctions 

φη (2) and $y, (2) are defined in the model atmosphere but 
the physical meaning of contributions from those terms is 
unclear. The contributions are included in the first two 
rows and the last row of each matrix. Second, only portions 
of eigenfunctions $5 (2) and φμη (2) are defined in the model 
atmosphere. To describe the incomplete sides of both eigen- 
functions an assumption must be made about the value of > 
at the surface and One] at the top of the atmosphere. Last, 
the equations derived for the general elements of each matrix 
are more complex than the formulas for the model FEM-C. The 
complexity of these equations makes it much more difficult 
to evaluate the elements in the first two and last two rows 
of each matrix. 

Assumptions are made and procedures are developed in 
an attempt to resolve these problems. In this model the mean 
state variables, u and T, are defined only at the N staggered 
levels. However, u and T values defined at the nodal points 


of $1 (2) and 6$ (2) are important in the Galerkin formulation 


N+2 
of the du/dZ and 3T/3Z terms. In my thesis experiments, the 
values of u and T are defined at the surface and top of the 

atmosphere; they are not defined at the nodal points of $4 (2) 


and 6$ (2) becuase these points are outside the model atmos- 


N+2 
phere. For constant shear with height, u and T are defined 
at the boundaries such that the shear in the two half layers 


at the boundaries is the same as the shear in the other layers. 
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The contributions from the perturbation quantities defined 

at the nodal points of $1 (Z) and ny C2) are not included in 
these experiments. To be consistent with the other N-layer 
models, all perturbation quantities are defined as column 
vectors of length N+l even though all quantities, except ver- 
tical veloc$to are defined at only N levels. The values of 
each perturbation at its lowest staggered level are stored 

as the second element in its column vector; zero is stored 

as the first element, which corresponds to the contribution 
from φη (2). To evaluate the staggered eigenfunctions defined 
and Z 


in the layers between the surface and Z and the 


2 


top of the atmosphere, it is assumed that the value at the 


Nara 


boundaries of those eigenfunctions is one-half. Thus, three- 
fourths of the eigenfunctions $^ (2) and Ong C2) are defined 
in the model atmosphere. 

The equations for the general elements of the five 
matrices are evaluated by substituting into equations (2-56)- 
(2-57) and (2-72)-(2-74) the formulas for $:,102, $; (Z), 
$: 4 (02, ΙΣΩΣ Y; (Z), ψ..1 (2), and VU: 2142) terms of the 
local coordinate Ë = Z - The equations for these eigen- 
functions defined for the levels 1, 2, i, and N-*1 are listed 
in Appendix E. The equations for the elements of the first 
two rows and the last row of the matrices P, 9, and R are 
derived using the formulas for levels 1, 2, and N-*l. However, 
the corresponding rows of matrices M and N were noro m 


lated using these formulas because the model results for the 
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barotropic experiment were not constant with height when the 
elements in those three rows of both M and N were calculated 
using these formulas. The reason for this is not clear. The 
equations for the general elements of M and N are simple 

enough to allow the terms in the affected three rows of each 
matrix to be deduced by reason. The equations for the elements 
of all five matrices are listed in Appendix F. 

The prognostic matrix equations for vorticity, diver- 
gence and temperature, listed in Appendix D, are simplified 
in a manner similar to the method described for the FEM-C 
model. The final form of the forecast equations for the 
vorticity and divergence vectors are the same as for model 
FEM-C, equations (2-66)-(2-69). The thermodynamic vector 


equations are 

















* E 
Tl |, ^ TlBAR, , * 2At.(Ql- u-M ^-N-T2 
* ο πα ου μι πα]... (2-75) 
— ΕΠ 
= ° e ul. ο 
T2 |, - T2BAR, , * 2at:[Q2 «iw M 1-N- T1 
“ἢ 
ο TE M s s; n "-R-W2] . (2-76) 
-------- —— n 


The six forecast equations are solved at levels two to ΝΑΙ. 
The surface geopotential and the diagnostic variables are 


calcuiated using the corresponding equations in model FEM-A. 
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3. FEM-B 

The FEM-B model defines vertical velocity, potential 
temperature, mean state potential temperature and diabatic 
heating at the unstaggered levels in terms of the eigenfunc- 
tions ο. The other variables are defined at the staggered 
levels in terms of the eigenfunctions $4 (2). The eigenfunc- 
tions are the same as defined for the FEM-A model, shown in 
Figure 3. 

The finite element approximations for the vorticity, 
divergence and thermodynamic equations are.derived by substi- 
tuting the eigenfunction expansion for each dependent variable 
into equations (2-35)-(2-40). The vorticity and divergence 
equations are multiplied by $; (Z) and integrated with respect 
to Z from the bottom to the top of the atmosphere. The resul- 
tant Galerkin formulation of the vorticity and divergence 
equations are the same as those derived for model FEM-A. The 
matrices in those equations, M, N, and P, are the same as 
defined for FEM-A, equations (2-56)-(2-57) and (2-73) NN 
thermodynamic equations are multiplied by V; (2) because poten- 
tial temperature is defined at the unstaggered levels. As 
before, the equations are integrated through the depth of the 
atmosphere. The resultant equations are listed in Appendix 
G. Four additional matrices are defined for the two thermo- 


dynamic equations. The mass matrix S is 


SM I / V4 (2) V, (2) a2 (2-77) 
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du. 


RucEmatrix O 1S defined for terms of the form az v 
T 2 ito . ^ 
g TEM Uu. $.(2)0,(2)y, (2)a82 . (2° 708) 
= jei-1 kei-1 Jz, ΤΘ; 
| | oT 
The matrix R is defined for terms of the form -7-W, 
Z 
TE a Ι'γτσν Ἱ οὖν 
R = ) ) T. f s= V. (2) yi (2) 82 à (2-79) 


OE Amp 


i 64 (22 (2) 5 (2) dz i (2-80) 


l| 3 

ii 
--η 
LJ 


As discussed in the FEM-A model description, the 
Staggered finite elements present problems for evaluating the 
elements of the matrices. In this model u is defined at the 
surface, the top of the atmosphere, and at the N staggered 
levels. The mean state temperature, T, is defined at the 
unstaggered levels so special definitions for it are not needed. 
Also, the contributions from the perturbation quantities 


defined at the nodal points of >, (4) and 6 (Z) are not in- 


N+2 
cluded in the model FEM-B model experiments. The staggered 
eigenfunctions, $; (Z), are evaluated at the boundaries using 


the assumptions discussed in the previous section. 
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The elements of matrices S, Q, R, and T are evaluated 


by substituting formulas for 6$. 


ο ου ο ο ο ο 


πος ψι (2) and Vi 102, defined in terms of the local 
coordinate n - Z -Zis into equations (2-77)-(2-80). Formulas 
for these eigenfunctions are listed in Appendix H. The equa- 
tions for the matrix elements are listed in Appendix I. 

‘The forecast matrix equations for vorticity, divergence 
and temperature are simplified in a manner similar to the 
method described for model FEM-C. The final form of the vor- 
ticity and divergence vector equations are the same as for 
model FEM-C, equations (2-66)-(2-69). The thermodynamic vec- 


tor equations are 

















τ τ. 
D = TlBAR, | t 2At[OI -p «Sg TET 
+ (Z) +57 *-Q-V1 - s l.np.w1] š (2-81) 
------- - ΟΕ 
i = IL 
dU T2BAR 1 + 2At{Q2 +u°S "T: Tl 
n s lo.v2 "osse W (2-82) 
C— maa mes ë ë Ππ| nD 











The vorticity and divergence equations are solved at levels 
two to N+1. The thermodynamic equations are solved at all 
unstaggered levels. The surface geopotential and the diagnos- 
tic variables are calculated using the corresponding equations 


in model FDM-B. 
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III. EXPERIMENTS AND RESULTS 


Three experiments are performed with each model; an 
initial perturbation in the meridional flow, flow over moun- 
tain topography and flow with a diabatic heat source. The 
first two experiments are performed with each model defined 
with six and then sixty layers. The heating experiment is 
repeated with six, twelve and sixty-layer models. The analytic 
solution of each experiment has not been derived. For each 
experiment, the 60-layer model results are intercompared to 
determine if the models are converging to the same solution. 
The standard of comparison for each six-layer model is its 
corresponding 60-layer solution. Temperature and divergence 
profiles are examined in each experiment. 

Several parameters are defined identically in each experi- 
ment. The vertical coordinate, Z, 1s defined between zero 
and one (1000-368 mb) and the vertical levels are equally 
spaced. The x-wavelength is 4,000 kilometers. The time step 
1s 17.7 minutes. The Coriolis parameter is defined at 45 
Eees latitude. There is no vertical shear in the u field 
and u = 10.0 meters/second (m/s). The mean state potential 


temperature increases with height from its surface value of 


EU Kelvin (K). 


A. ROSSBY WAVE EXPERIMENT 
Rossby waves are generated in each model using an initial 


perturbation, v’' = 5.0 m/s, in the cosine term of the 
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meridional flow. All other perturbations are initially zero. 
There is no diabatic heat source and no mountain topography. 
The latitudinal variation of the Coriolis parameter, 8, is 
defined at 45 degrees latitude. The forecast length of this 
experiment is 96 hours. 
l. Sixty-Layer Models 

"The 60-1layer FDM-A, FDM-B, FDM-C, and FEM-C models 
converge to the same temperature and divergence solutions 
(Figures 4-7). It should be noted that the phase of each 
variable is defined between zero and 360 degrees. There is 
a discontinuity in the phase profile if the phase passes 
through zero degrees. The height at which the temperature 
phase discontinuity in model FDM-A occurs differs from the 
other three models because temperature is defined at the 
staggered levels in FDM-A. The four models represent the 
same physical solution, which is called the consensus solution. 

The FEM-B temperature amplitude is slightly larger 
than the consensus amplitude (5%) and an amplitude oscillation 
1s present in the lowest three layers of the atmosphere (Figure 
8). The jagged profile may be caused by the terms in the 
matrices which represent contributions from the eigenfunctions 
near the lower boundary of the model. There is a discontinuity 
between the temperature phases in the lowest two layers of the 
FEM-B profile (Figure 9). Apart from this feature, the 
general shape of the phase profile is similar to the consensus, 


although the phase is five degrees less and passes through 
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Figure 4. Sixty-layer Rossby-wave experiment at 96 
hours. Temperature amplitude profiles are 
compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C: 
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TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 
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Figure 5. As in Fig. 4 but for temperatures 
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DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 
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As in Fig. 4 but for divergence amplitude. 
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Figure 7. As in Fig. 4 but for divergencegpassec 
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Figure 8. Sixty-layer Rossby wave experiment at 96 
hours.  FEM-B temperature amplitude profile 
is compared with the temperature amplitude 
profile of FDM-C, which represents the 
consensus profile. 
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TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 
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Figure 9. As in Fig. 8 but for temperature phase. 
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zero at a lower level than the consensus. The divergence 
amplitude of model FEM-B is within five percent of the con- 
sensus (Figure 10) but its phase profile (Figure 11) is quite 
different from the consensus. 

The model FEM-A temperature and divergence ampli- 
tudes vary significantly from the consensus. The temperature 
amplitude profile (Figure 12) has a large amplitude oscilla- 
tion in the lowest three levels above the surface, possibly 
caused by the choice of eigenfunction representations used 
to evaluate elements of the first two rows of the matrices. 
The FEM-A temperature amplitude is larger than the consensus 
throughout the atmosphere. The temperature phase passes 
through zero at a lower level and is five degrees less than 
the consensus (Figure 13). The maximum FEM-A divergence 
amplitude occurs near the middle of the atmosphere and its 
magnitude is ten percent less than the maximum of the consensus 
(Figure 14). The consensus divergence phase is nearly con- 
Stant with height, but the phase profile for model FEM-A 
steadily decreases with height (Figure 15). 

e Six Layer Models 

The comparison of six and sixty-layer profiles for 
variables defined at staggered levels may be initially mis- 
leading. The first staggered level in a six-layer model occurs 
at Z = 0.0833. The lowest staggered level ina 60-layer 
model is defined at Z = 0.0083, which may be mistaken for the 


surface in the graphs. When the values of a six-layer model 
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Figure 10. As in Fig. 8 but for divergence amplitude. 
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Eure 11. As in Fig. 8 but for divergence phase. 
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Figure 12. As in Fig. 8 but for model FEM-A compared 
with the consensus, FDM-C. 
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TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 
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Figure 13. As in Fig. 12 but for temperature phase. 
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DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 
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Figure 14. As in Fig. 12 but for divergence amplitude. 
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DIVERGENCE PHASE FOR V = 5.0 CASE (60 LAYERS) 
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ENcure 15. AS in Fig. 12 but for divergence phase. 
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coincide with a 60-layer profile, the models are considered 
to represent the same physical solution, even though the 
six-layer model has a smaller vertical domain for staggered 
variables. 
The temperature amplitude profiles of the six-layer 

FDM-A, FDM-B, FDM-C and FEM-C models are identical with their 
corresponding 60-layer results (Figures 16-19). The previously 
discussed problems in the lowest layers of models FEM-A and 
FEM-B are quite evident in their six-layer profiles. The 
temperature amplitude of FEM-B is virtually identical with 
its 60-layer profile above the first two layers (Figure 20). 
The six-layer FEM-A profile (Figure 21) does not agree with 
its 60-layer profile below Z = 0.50 (600 mb). The large 
amplitude oscillation in this profile destroys the integrity 
of the solution in a significant portion of the atmosphere. 

| The six-layer divergence amplitude profiles for the 
grid C models are identical with each other (Figures 22-23), 
and the divergence profiles for the grid B models are nearly 
the same (Figures 24-25). All four models approximate their 
corresponding 60-layer solutions with similar accuracy. The 
Six-layer FDM-A model does not approximate the curvature of 
the consensus 60-layer divergence profile as well as the 
grid B and C models (Figure 26). The six-layer FEM-A model 
has difficulty approximating the location and magnitude OE 
the mid-atmosphere divergence maximum of its 60-layer solution 


(Figure 27). 
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Figure 16.  Six-layer Rossby wave experiment at 96 
hours. Temperature amplitude profiles are 
compared for the six-layer and 60-layer 
FDM-A models. 
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Figure 17. As in Fig. 16 but for model FDM-B. 


56 





— À PR e T rm PPM 


— —— —— . 


0.9 


0.8 


0.7 


Q.6 


0.5 


HEIGHT, Z = — LN(P/PO) 


0.5 


Q.2 


0.00 0.02 0.04 0.06 0.03 0.10 0:12 0.14 


TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FDM-C) 


TEMPERATURE AMPLITUDE (KELVIN) 


LEGEND 
- FDM-C 6 LAYERS 
^ — FDM-C 60 LAYERS 


390 
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τε 15. AS in Fig. 16 but for model FDM-C. 
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Figure 20. As in Fig. 16 but for model FEM-B. 
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Figure 21. As in Fig. 16 but for model ere a 
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Pagure 22. As in Fig. 16 but divergence amplitude of 
ΕΡΕ. 
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Figure 23. As in Fig. 22 but ftorumode ΗΕ CF 
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Figure 25. As in Fig: 22 bu o model FEM-B. 
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Figure 26. As in Fig. 22 but for model FDM-A. 
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Figure 27. AS in Fig. 22 but for medel πρ 
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B. MOUNTAIN TOPOGRAPHY EXPERIMENT 

The forced vertical velocity term, MTS, in the surface 
geopotential forecast equation (2-3) is non-zero in this 
experiment. It represents the contribution to surface geo- 
potential from air flowing over mountain topography which 
varies sinusoidally in the x-direction and has no variation 
in the y-direction. The mountain ridge-to-valley height 
difference is 1,500 meters. To reduce the trauma for the 
model, the mountains are gradually "built" to their full 
height over a period of 36 hours. Thus, the forced vertical 
IE τν increases in the first 36 hours of the forecast 
period and is constant for the remainder of the 96-hour fore- 
cast period. The equations used to define the forced verti- 
cal velocity are included in Appendix J. Beta and all initial 
perturbations are zero in this experiment. 

1.  Sixty-Layer Models 

The FDM-A, FDM-B, FDM-C, and FEM-C models converge 
to the same physical solution. The amplitude and phase pro- 
files of temperature and divergence are depicted in Figures 
mo-31. 

The FEM-B model again has a jagged temperature ampli- 
tude profile in the lowest two layers (Figure 32). The 
temperature amplitudes in the remainder of the atmosphere are 
within 5% of the consensus. The FEM-B temperature phase has 
oscillations in the bottom three layers and the top two layers 


of the profile (Figure 33), but the rest of the profile is 
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Figure 28. JSixtyslayer mountain topography experiment 


at 96 hours. Temperature amplitude profiles 
are compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C. 
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Figure 30. As in Fig. 28 but for divergence ampuiss e 


70 


DIVERGENCE PHASE FOR MOUNTAIN CASE (60 LAYERS) 


0.9 


HEIGHT, Z = — LN(P/PO) 


0:5 


ox 


—— 141 


0-17 


0 50 100 150 200 250 300 350 400 
DIVERGENCE PHASE (0-360 DEGREES) 


LEGEND 
azFDM-A ^' 
ç = FDM—B 
4 2 FDM-C 
+ = FEM—C 


]Eicure 31. As in Fig. 28 but for divergence phase. 
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Figure 32.  Sixty-layer mountain topography experiment 
at 96 hours.  FEM-B temperature amplitude 
profile is compared with the temperature 
amplitude profile of FDM-C, which represents 
the consensus profile. 
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within one degree of the consensus phase. The FEM-B 
divergence amplitude and phase are identical with the 
consensus values. 

An oscillation also exists in the lowest three layers 
of the temperature amplitude profile for model FEM-A (Figure 
34). The amplitude is 30$ higher than the consensus near the 
surface and this difference decreases with height. The tem- 
perature οσα are nearly identical. The FEM-A divergence 
amplitude profile is similar to the consensus, but the amplitude 
does not decrease as fast with height as the consensus 
(Figure 35). The divergence phase profile for this model is 
the same as the consensus. 

2s. Exe Bayer Models 

The temperature amplitude profiles of the six-layer 
models FDM-A, FDM-B, FDM-C and FEM-C are identical with each 
other and also with the consensus solution (Figures 36-39). 
The six-layer staggered finite element models are nearly 
identical with their 60-layer solutions above the first two 
layers (Figures 40-41). The previously identified inadequate 
representations of the temperature amplitude in the lowest 
layers exist in this experiment. All six models approximate 
their 60-layer divergence amplitude profiles quite well 
(Figures 42-46), however the FEM-A model (Figure 47) does 
not approximate the curvature in the lower portion of its 


60-layer profile as well as the other models. 
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Figure 34. As in Fig. 32 but for model FEM-A compared 
with FDM-C. 
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Six-layer mountain topography experiment at 
96 hours. Temperature amplitude profiles 
are compared for the six-layer and 60-layer 
FDM-A models. 
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36 but for model FDM-B. 
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Figure 38. As in Fig. 36 but for model FDM-C. 
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Figure 39. As in Fig. 36 but for model FEM-C. 
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Figure 40. As in Fig. 36 but for model FEM-A. 
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Figure 41. As in Fig. 36 but for model FEM-B. 
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τε 42, AS in Fig. 36 but for divergence amplitude. 
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Figure 43. As in Fig. 36 but qi I u amplitude 
for model FDM-B. 
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Figure 44. As in Fig. 36 but divergence amplitude 
for model FDM=C. 
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Figure 45. AS in Fig. 5 Dur EE divergence amplitude 
for model FEM-B. 
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Figure 46. As in Fig. 36 but for divergence amplitude 
for model FEM-C. 
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Figure 47. As in Fig. 36 but for divergence amplitude 
for model FEM-A. 
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E; DIABATIC HEATING EXPERIMENT 
A diabatic heat source is defined in the layer between 
Z = 0.40 and Z = 0.60 (670-549 mb). The rate of heating is 


constant.in time and varies in x and Z, 


(Z-Z) 


Q(x,Z,t) -= HEATING:cos [z te 


— T] cos(ux) , -ῃ 
U 


L) 


v 


where HEATING is 5.0 K/day, Zu is the midpoint of the heated 


layer, and Σι and ZU are the lower and upper boundaries, 


respectively, of the heated layer. The diabatic heating 


--,. 01 and 02, are defined in the initialization subrou- 


tine and stored for use in the forecast subroutine. Beta and 
all initial perturbations are zero. The forecast length is 
ime hours. 


im  olxtyrlayer Models 


For the diabatic Meat ine function defined in equation 
(3-1), the maximum heating occurs at Z = 0.50, the midpoint 
of the heated layer. The models defined using grids B and C 
define temperature and the heating functions at this point. 
The grid A models do not have temperature and diabatic heating 
defined at this point so the maximum rate of heating in these 
models is slightly less than in the other models, and the 
maximum heating occurs throughout one layer rather than 
occurring at one point. The heating rate at each level is 
listed in Appendix K for the six, twelve and sixty-layer 


staggered and unstaggered models. 
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The 60-layer profiles for the six models are quite 
Similar, but the differences occur because the models are 
responding to different forcing. The temperature amplitude 
profiles-for the B and C grids come to a sharp point at 
Z = 0.50 and the grid A models have a square-nosed profile 
around this point (Figures 48-50). The height in the atmos- 
phere at which the model experiences the onset of forcing is 
slightly different between the grids. This explains the small- 
scale oscillation in the temperature amplitude profiles near 
the boundaries of the heated layer, Z = 0.40 and Z = 0.60. 

The previously identified temperature amplitude oscillations 

in the lowest layers of models FEM-A and FEM-B are not evident 
in this experiment. The temperature phase profiles are nearly 
identical, except that the grid A models do not have απ ο 
of a spike at the boundaries of the heated layer as the other 
grids (Figures 51-53). In summary, the 60-layer temperature 
profiles of all six models represent the same physical response 
to the diabatic heating. 

Divergence is defined at the midpoint of the heated 
layer only in the grid C models. Consequently, the shape of 
the divergence amplitude profile for the grid C models is 
different than the other four models. The grid C profiles 
have a sharp point at 2 - 0.50 (Figure 54). However, the 
minimum divergence is not symmetric around Z = 0.50 in the 
grid A and B models (Figures 55-57). The divergence ampli- 


tude is identical outside the heated layer for all models 
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Figure 48. Sixty-layer diabatic heating experiment at 
12 hours. Temperature amplitude profiles 
are compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C. 
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Figure 49. Sixty-layer diabatic heating experiment at 
12 hours. FEM-A temperature amplitude profile 
is compared with the temperature amplitude 
profile of FDM-C, which represents the 
consensus profile. 
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Figure 50. As in Fig. 49 but for model FEM-B compared 
with the consensus, FDM-C. 
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48 but for temperature phase. 
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Figure 54. Sixty-layer diabatic heating experiment at 
12 hours. Temperature amplitude models are 
compared for the two grid C models. 
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Figure 55.  Sixty-layer diabatic heating experiments 


12 hours. Temperature amplitude profiles 
are compared for models FDM-A and FDM-B. 
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Figure 57. As in Fig. 55 but models FEMA ana ee 
are compared. 
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except FEM-A (Figure 57). The divergence phase profiles are 
Erruatlvoydentical for models FDM-A, FDM-B, FDM-C and FEM-C 
(Figure 58). The divergence phase passes through zero at a 
different level for model FEM-B than the consensus, but this 
difference is not physically significant (Figure 59). The 
FEM-A profile differs slightly from the consensus outside the 
heated layer (Figure 60). In conclusion, the six models are 
converging to the same 60-layer divergence solution inside 
the heated layer, but model FEM-A does not agree with the 
solution of the other five models outside this region. 
2. Six and Twelve-Layer Models 

The difference between grids is more evident in this 
experiment than in the other experiments. Each model is run 
with both six and twelve layers and the results are compared 
with the 60-layer consensus of the six models. The tempera- 
ture amplitude consensus profile has a small-scale decrease 
at each boundary of the heated layer that cannot be reproduced 
using either six- or twelve-layer resolution. 

The six-layer grid A model temperature and divergence 
fields barely respond to the diabatic heating. The grid A 
models have identical temperature responses for both six and 
twelve layers. The six-layer perturbation, constant with 
height in the heated layer, is an order of magnitude smaller 
than maximum consensus perturbation (Figures 61-62). The 
grid A models have a stronger response to the heating with 


twelve layers than with six layers, but the maximum amplitude 
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Figure 58. As in Fig. 48 but for divergencegbbsscs 
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Figure 60. As in Fig. 57 but for divergencerphases 
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Six-layer diabatic heating experiment at 12 
hours. Temperature amplitude profiles are 
compared for the six-layer and 60-layer 
FDM-A models. 
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is only 60% of the 60-layer consensus response (Figures 
63-64). For both the six- and twelve-layer profiles, the 
temperature amplitude at the two staggered nodal points 
defined closest to the heated layer lie directly on the 60- 
layer profile. This may indicate that for much less than 
60-layer resolution the amplitude of the spike will not equal 
the maximum amplitude in the 60-layer experiment. The 
divergence fields of the six-layer grid A models is also an 
order of magnitude lower than their corresponding 60-layer 
solutions (Figures 65-66). The characteristic sharp divergence 
decrease in the heated layer 1S roughly approximated by the 
12-layer models (Figures 67-68). The minimum divergence is 
much larger than in the consensus profile and the base of the 
Spike decrease 1s much broader than in the 60-layer solutions. 
Outside the heated layer, the shape of the 12-layer profile 
1s πίτας to the consensus, but the amplitude is lower. 

The six-layer temperature amplitude response of the 
B and C grid models are virtually identical (Figures 69-72). 
The base of the profile spike is broader in the six-layer 
models than in the consensus, but the four models closely 
approximate the magnitude and width of the tip of the spike. 
The 12-layer temperature profiles for grids B and C are also 
identical (Figures 73-76). The profiles slightly underesti- 
mate the magnitude of the perturbation in the top half of the 
Spikes, but the width at the base of the spikes are the same 


as in the 60-layer consensus. The small scale dip at the 


107 


TEMPERATURE AMPLITUDE FOR HEATING CASE (FOM—A) 


0.9 


0.8 


HEIGHT, Z = — LN(P/P0) 
O 





0.3 
0.2 
"m 
l 
κ 
ο cC" cn E E 
9 0.25 | 0.50 0.75 i i25 1.59 1.75 2 dis 


TEMPERATURE AMPLITUDE (KELVIN) 


LEGEND 
C 2 FDM-A 60 LAYERS 
c = FDM-A 12 LAYERS 


Figure 63.  Twelve-layer diabatic heating experiment at 
12 hours. Temperature amplitude profiles 
are compared for the 12-layer and 60-layer 
FDM-A models. 
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Figure 65. Αη" ΥΕ» divergence amplitude. 
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Figure 66. As in Fig. 62 but for divergence amplitude. 
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Figure 67. As in Fig. 63 but for divergence amplitude. 
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boundaries of the heated layer in the consensus profile is 
not reproduced by the 12-layer models, otherwise the profiles 
agree with the consensus outside the heated layer. 

The diabatic heating causes different divergence 
amplitude responses for the B and C grids. The grid C models 
have similar six-level profiles (Figures 77-78). The minimum 
divergence, which occurs at the center of the heated layer, is 
nearly the same in the six and Sixty-layer models. There are 
differences between the two grid C models away from the diver- 
gence minimum. Outside the heated layer the six-layer 
divergence amplitude of the FEM-C model provides the closer 
approximation of the 60-layer consensus profile. Similarly, 
the 12-layer FEM-C profile (Figure 79) is a better approxima- 
tion of the consensus than the 12-layer FDM-C model (Figure 
0). 

Fór both six and S layers, the grid B models very 
poorly approximate the 60-layer divergence amplitude within 
the heated layer. The 12-layer profiles (Figures 81-82) 
are identical for the grid B models and are closer approxima- 
tions of the consensus 60-layer profile than the six-layer 
results. The 12-layer model is asymmetric within the heated 
layer and does not have the dramatic decrease in divergence 
that exists in the consensus. The twelve and sixty-layer 


profiles are similar outside the heated layer. 
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Figure 77. As in Fig. 61 but divergence amplitude for 
model FDM-C. 
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Figure 79. As in Fig. 63 but divergence amplitude 
for model FEM-C. 
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Figure 82. As in Fig. 63 but divergence amplitude for 
model FEM-B. 
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IV. CONCLUSIONS 


Although the analytic solutions of each experiment have 
not been derived, some conclusions about the characteristics 
of each model can be drawn. During the Rossby wave and moun- 
tain topography experiments the staggered finite element 
Models, FEM-A and FEM-B, display unusual temperature ampli- 
tude behavior in the lowest two layers of the model. The 
oscillation in the temperature profiles of both models may 
be generated by the matrix elements which represent the 
contributions from the eigenfunctions near the surface. It 
is possible that better representations of the staggered 
eigenfunctions near the lower boundary may reduce or eliminate 
this problem. Jagged eer seus profiles were not observed 
in the diabatic heating experiment. The 12-hour forecast 
period in the heating experiment may be too short for the 
profile discontinuities to grow to substantial amplitudes. 

The differences between the grids is most apparent in the 
diabatic heating experiment. The 60-layer temperature pro- 
files of all six models represent the same physical solution. 
The FEM-A model has a slightly different 60-layer divergence 
response outside the heated layer than the other five models. 
The grid B and C models have identical temperature amplitude 
profiles for both six and twelve layers. However, the diver- 


gence amplitude profiles are quite different and the grid C 


E29 


models provide the better approximation for both resolutions. 
The grid A models barely respond to the diabatic heating with 
Six layers, and poorly approximate their consensus profiles 
with twelve layers. There are indications that both grid A 
models may not converge to their 60-layer temperature and 
divergence solutions with much less than 60-layer resolution. 

In all experiments, the grid B model produces viru u 
identical results, aside from the small amplitude oscillation 
in the FEM-B temperature profiles. Based on these experiments, 
no accuracy is gained by using the finite element approxima- 
tions with grid B. In two of the three experiments the grid 
C results are identical. The FEM-C model provides a closer 
approximation of the consensus 60-layer divergence profile 
than the FDM-C model in the diabatic heating experiment. 
Generally, the lower resolution FEM-A model does not approxi- 
mate its 60-layer solution as well as the FDM-A model. It is 
not known if the different approximation characteristics of 
the FEM-A model are due solely to its profile problems at the 
bottom of the atmosphere. 


The importance of the φη (2) and φ (Z) eigenfunctions in 


N+2 
the staggered finite element models is not clear. They are 
prominent terms in the mathematical derivations, but their 


physical significance is not apparent. 
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APPENDIX A 


FINITE DIFFERENCE APPROXIMATIONS 


1. For terms of the form S W: 
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APPENDIX B 


GALERKIN FORM OF FEM-C PROGNOSTIC EQUATIONS 
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APPENDIX C 


MATRIX ELEMENTS FOR FEM-C 


` . — - = «πε ' 
1. Notation: Δ. δι 21.1, but Δ» Z, - 2; 
! = ! σι 
A 2 z m 
= = -- = ΕΠ : 
πο. ο. = u, T at 1 vertical level 
1 B 1 
2. Matrix decrnit ion; M (row, column) 


EI 
The 1 row of a matrix contains the weighting factors of 
the terms that affect the value of the forecast variable at 
th : ! : ; th 
the 1 vertical level. The weighting factor in the j 
column of the 22 row is the influence on the forecast varia- 
En : Γη I 
ble at the 1 level from the variable at the j vertical 


level. Level i = 1 is the surface, level i = N+l is the top 


of the atmosphere. 


3. To calculate the terms of each matrix, substitute in the 

matrix formula For φ. μη (2), o, (2), Φα (2) using functions 

defined in local coordinates with respect to >; (2). Define 
= _7! 

£ TEN Ζ Zs 


The eigenfunctions in local coordinates are 
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matrix elements corresponding to the surface 


and the top of the atmosphere, use the above formulas for 


i = 1 and i = N-*1, respectively. 
Mass matrix, M: 
a. General formulas, 
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C. At level i = N+t+l: 
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APPENDIX D 


GALERKIN FORM OF FEM-A PROGNOSTIC EQUATIONS 


1. The vorticity equations, (2-35)-(2-36), have the same 
form as the vorticity equations in FEM-C. While the FEM-A 
equations look the same as equations (B-1)-(B-2) in Appendix 
B, the equations for das $. 


J 
$(2) eigenfunctions are defined for the staggered levels in 


and $4. are different because the 


model FEM-A. 


2. Divergence Equations (2-37)-(2-38): 
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APPENDIX E 


LOCAL COORDINATES FOR STAGGERED ó; (2 


1. To evaluate the integrals in the finite element approxima- 
tions of the vorticity and divergence equations in models 


FEM-A and FEM-B, express $; ), δι» Φι.1» Vi, Vio Vul] 


and Ui; 5 in terms of the local coordinate £ = Z2 -2;. The 


general form of the local coordinate system is shown below. 


n 
e e ae ee e or 
; AT 
Zi $75 Hs 
D Eo Φ 
; EE Y 
eii SES ee 
κ EUER p. 
ιο 
De In terms of local coordinates fos 0, (5), the equations 


tocit Pr ee ee and ψᾳ.2 are listed below 


The notation used is defined in Appendix C. 
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See With respect to the function >, (8), the layer of interest 


1s Zo < 2 < 2 


jM 
Or 282 


eres. 5) 
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Eigenfunctions $, (&), 


o5(8), Wy (8) and Wo (8) are defined in this layer. 
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Eigenfunctions dy: 95 $4! 
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55 
Ë + A? N 
2 2 
AS paces D 
D E) ΚΞ 
2 -E +A, 
-- ο. 
: -T 
(E) Es 
Y tE M eq 
j A, 
1 
aS BEN a 
ve). = 
E _ š HIA) A. A? 
πο ο 
Δ} 2 — — 3 
JL 
net A! 
_ 2 2 2 
ET pu a. 


pee With respect to the function 94 "tne layer of 


NEE 


9 . I . * 
' E AG 
interest is Zu < 2 < Zed Or ANGI NONE 25N4l' |Ergenrtunctrons 


ON4-2 ' φμει’ ὃν, Pyare YN and ψμ-1 3τε defined in this layer. 


147 


[> 


N+1 


-&6 4 t 
^ Nl 
Lt 


NES 


dors 


SY 


Naw 
; 
NFL 





N+1 


μη 


|^ 


|^ 


|^ 


|^ 


|^ 


[^ 


Pa 
WD yet 


Punel 


€ 9 
z 


z 
' 
aes 


ÔN+1 





Τῇ 
2^N*1 


JYI 
| ^ 
C 


JYI 
|^ 
O 


Y 
|^ 
e 





APPENDEAX F 


MATRIX ELEMENTS FOR FEM-A 


e Mass matrix, M, using the notation defined in Appendix C: 


a. General formulas, i = 2,...,N: 
M(1,i1-1) = τ (F-1) 
orc]: ΕΠ 
"NIE MET (F-2) 
1 3 151] B 
PEE =A (F-3) 
; 6 i-1 


b. At level i = 1: 


ye ME 
M (15 D) B 3^2 (F-4) 
MPO NEC. (F-5) 
á |» 6 2 f 
c. At level i = N+1l 
M(N+1,N) = TA (F-6) 
: EE SON EI 
ZU NT EU MES ο} D 
᾽ 3 ΝΑΙ 2 N+1 
2. Matrix N: 
a. General formulas, i = 2,...,N 
Ri eye a aa) (F-8) 
4 jo T i-l i 
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N(i,1) = 12 Cy EU Δι) 
+ =u m +A.) F-9 
{αι δεν] δι (5-5) 
Niir s OH CUOI M (F-10) 
: 1255938] X i+l 
b. At level i= 1: 
f FERT l,- | ' 
N(1,1) = y5(¥74,) * z(u,- A.) (Ε-1ΙΝ 
ο. πα πα... (F-12) 
το ανα» 
C. At level 1 = N+tl: 
η, πο ο πο ο... (F-13) 
, 12 Nl N 
Na =a S a a MP ME 
Ç 24 ` N+2 Ντι 12 N+l N+l 
e ae ] 
μι C sn πει) (FD 


s Ma ri È: 


a. General formulas, i -» 3,...,N 
Por 3-2) = (αν ταις 1) le Δι 1 
(a1)? 
= — πα (F-15) 
° S 1 
48 (Δι) (A; 1) 
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BP(i,i-l) = {(u. 








DT au Ai I A; Ai 
go = ate. πμ ορ - b 
1 iel 2 O° As 1 AEG B+ A; 
(at)? Ue (aij 
=F s ae + ma sms l aros 
48: (Δ:) E 
y A (at)? 
ory) = itu. -u._,)-(>= 3 
i i-l Δ. 4g (A,)* 
3« A! Eye 
αρ τας) ga ` p 
L 1 E 
itl 48: (A +1) 
is sta E it 
"μ.ο... 
To DUC (At)? 
1 itl 7 
se) 2 7^ 2 P 
1151 itl 8- (A, 1) 48- (Δ. τη) (Ai) 
(F-17) 
2 
A. A! (A!) 
E. - _ + 
itl itl 8*A: 1175141 
(at)? 
--------------- (F To) 
JL MD 
b. At level 1 5 1: 
P(l,l) = $(ü5 uj) (F-19) 
EINEN ig(u, -u) (F-20) 


5 


c. At level i = 2: 











' i 2 
12 9008 nS as (as) 
ον i) = 48 (u, - u) + (u, - u>) : [SA Pe — Έτους 
3 48+ (AQ) 
A A! A! 
PEN A η... M NEN 
3 3 3 
2 
(As) (as)? 4” 
4- 5 - + 5 ] (F-2 
48 (A4) 8:A,*A3 — 48(A4) (A3) 
_ AB, aj as)? (as)? 
P (235) = uou sal ' ala . C CR M 
3 2 6.43 4-4: — 8:54:53 48 (54) ^ (a3) 
(F-23) 
C. At level i = Ntl: 
3 
u s: ee) 
μ.ο... Riel (F-24) 
Ntl CN /6-Ó& ao A 
ΠΤΙ Ν 
ΒΡΙΝΕΙΙΝ = (Ge, te.) EMI ΚΣ 
(N+1,N) = tag Yyy2 "ἄμμι Uyl Tn 2 7 Gear 
A: A. ο... 
QNM _ N+1 N+1 
QAM. TESTEN | 2 
N N-41 48 (Ava? 
2 3 
(ATE) (ΑΙ, 
= — a: — l (F-25) 
(Αι) (At) θα μοι 
M - 2 = AN+1 
P(N+1,N+1) = {38(0y+2 TYN) * (ἄμμι ταν) ἱπτη--- 
N41 
2 
[Απ 
- Ntl ἢ) (F-26) 
1850 d 
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4. Matrix Q: 


-wGeoneral formulas, 2 = 2,...,N 
ο πα αι ο) (F-27) 
i , e ys 
Ce) =e. “πα (F-28) 
, Bu sy 
(CERES l(u - u.) (F-29) 
Qii,l EN CUu 
b. At level i = 1: 
Q(1,1) = S -u.) ΠΠ 
, ο. αι 
(1.2) 5: πο (Ε-31) 
Qtt, "7 g5-U 
ee At level i = N+1: 
Ο(Ν41,Ν) mu (F-32) 
, EN COSS 
CILE το σα) (FS33) 


meee Matrix R: 
The formulas for the elements of this matrix are the same 
Ezmchose for matrix P, with m replacing SIE j r Ons 


i 15)-(F-26). 
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APPENDIX G 


GALERKIN FORM OF FEM-B PROGNOSTTCFEOUAT TONS 


l. The vorticity equations have the same form as the vor- 
ticity equations for model FEM-C, equations (B-1l)-(B-2). 


However, the equations for $4 7 Q. and φις are different 


7 
because $(Z2) is defined for the staggered levels in FEM-B. 


2. The divergence equations are the same as the divergence 


equations in model FEM-A, equations (D-1)-(D-2). 


3. Thermodynamic Equations (2-39)-(2-40): 


aT aj n 
2 get J 4020,02)82 
J=1-1 Z 
0 
τ μμ. - 
== u | ; ΤΠΕ b. (Z) (2) h o K 
jei-1 k-i-1 7 ^ E 
| | Z 
=) Y 2 uv (e) J L 6, (2) y, (2) aZ 
Ere ο — ! 
Rjsicl kei-1 2" 7, MEE 
2 
1151 UY - ak f T UE 
- ) TW (E) ——= ww, (2). (Z) dz 
eeu Zo dZ K 1 
i+l ! an 
+ Σ ott) f ps (2)y, (2) az Te 
--- i) T 
j= 1=1 Zo 
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111 ar? 2η 
j=i-1 dt z J PN 
0 
ic TEM ος on 
= TE E 2) ψ. 
f i+2 EPA AE am dé. 
, 0 
a d E on HE 
- oT W(t) ο ο σα 
j=i-1 k=i-1 7 4 A az RN Yi 
itl j on 
+ * EL ; J 94029, (2) dz (G-2) 
Τη 0 
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APPENDIX H 


LOCAL COORDINATES FOR UNSTAGGERED ψ. (2) 


l. To evaluate the integrals in finite element approximations 
of the thermodynamic equations in model FEM-B, express 9:427 
9:11" $i: $i 1, ψ 1’ Ws and Yil in terms of the local 


coordiante n = Z τς. The general form of the local 


coordinate system is shown below. 





i+2 
n = x 
i+1 i+1 i+1 
A. 
mc — ----. = —I+1 
Zia E 2 n 
2 n=0 τ 
BT 
z n = 7 Φ 
£i ? = - á, T 
q 
-1 
2. In terms of the local coordinates for v i (n), the equations 
for 9:12* 914,1" Vir Vilcl' pr ji and Vill are listed below. 


The notation used is defined in Appendix C. 
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m2 itl να αι 
τη ν΄... Ξε 
T2 
EUN | I 
$; (n) - A AT OS nos 54. 
D 
η ἘΞΔ' 
ET M. qu 
ES pU T 
$. (n) = 
ie ο τι. 
2 L ΠΠ Έλι 25 A! 
ΔΙ 2 i*l — 151 
1. 
n - =A. 
E 2 i*l s ' 
owa S eee ey 
n +A! 
V 0 = n b OP = 
1 
ψ.(π) = 
au 
τη. 
1-1 ' 
AT NEUE d 
απο A 
E n 
πη - DM Οι 


3. With respect to the function Y (n), the layer of interest 
is 24, <42< 2, or O<n< A5.  Eigenfunctions Pr 027 Q37 
y, and y, are defined in this layer. 
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Φ 
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Z — m = > 
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2500 oe n= =? 
2 2 e 
TTTTTTT ο 
Zo 0*0 Ha 
4 + 
D 
-n +54 
= ο E 
$4 (n) - a ο --σδ2 
ww 
τ» » .- 
As 0 < n < 555 
“2 
$5 n) - cue 
-7 +5 
2 2 lae A 
Δρ 2-2 = T = %2 
1; 
nm s ia) : 
DON EE À 252 SN 
7n * 4) 
p (n) = 5 ο ο... As 
ψρίη) = ra O<n< 4s 


4. With respect to the function y(n), the layer of interest 
15 45 5 2 « Za. Eigenfunctions $17 25! 94: $a Wye Yo and 


Ņ are defined in this layer. Formulas for Qar Ui wy, and 


V3 are obtained from the general formulas with i - 2. 
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5. With respect to the function Vey (ns the layer of interest 
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APPENDIX I 


MATRIX ELEMENTS FOR FEM-B 


ου. ος” and © are the same as those defined in 
model FEM-A. The equations for their elements are listed in 


Appendix F, equations (F-1)-(F-26). 


Kass matrix S is the same as matrix M in model FEM-C. 
EEEEcquatrons for the elements of S are the same as equa- 


tions (C-4)-(C-10) in Appendix C. 


bee Matrix R is the same aS matrix R in the model FEM-C. 


Equations (C-25)-(C-31) in Appendix C define the elements of 


IE: 


Z Matrix Q, using the notation defined in Appendix C: 


a. General formulas, i = 2,...,N: 
. C 
En 1-1) - πα. i m ol (I-1) 
48- (Δι) 
2 
S OD), UA ORAE) 
Q(i,i) - (u; -u, ):[———, + — 15 
24 (Δι) 16-(Δ.) 
2 2 
B E (E) ET) Ρα.) 
+ (u ο m + — 1+1 oi + 1+1 ] (I-2) 
in UU NATU EE ICE 48(A. .)2 
111 itl 111 


Tol 


(A!) J(AT 9m 
O (1;71+1) = (u u.) [m r 
νυ,’ 
111 itl 
2 2 
s | 
+ ee + ας μαι “itl 
48* (A; , i) 24- (A, ,) 
ha OD 
16- (Δ: γη) 
2 
CANCER 
ος " = Ti . Itl 
Ou 1356294 3 (ας νο τι μι) m d | ] 
ΠῚ 
Ὁ. At level i= 1l: 
πη 
" .. 
T (a5) * ú aa 
QD) eel t-bone ] 
UI mE ο ο (A) 
λ α. (0) * 
Q(1,3) = (u> ~- u3) | 
ΠΣ 
C. At level i = N+tl: 
2 
- _ (A 
Q(Nt1,N) - p p i aa t 
48* (A...) 
N41 


(P= s) 


(1-5. 


(1-9) 


(I-69 


(Ic 


(1339 


2 - 


13: (u -u,,,) ú _ CA Ὁ 
Q(N+1,N+1) = = ο 7 —7 Nee 2 
24: (À 
N+1 


(A ΠΟ 


16. (Δι, 1) 


gee Matrix T: 


v 


a. General formulas for 1 = 3,..: ,N 


E 2 Εν s 
IINE LE OECD 


E D — oA Uy ‘Toon * — 2474. 
L JE T 


' 2 $ f aa S 1 2 
O e AD o wur (ap -- 
132-41 T972 Δι γι 64* Ai 1 


EUN α upar ο) 
192-4; 1 64.åd, 48-A; 





E | m". 2 
: 1l- (Ai) . 7 (A η} (Δ1} N 1 7 (A; 1) | 
DIM js Ζ. 


2 i ' 
1) . 7 (A; 11) 043) 
1932*5; 1 48*45., 


' 2 ° 3 
17 (A;) TI (Ai, 


Πιο il 
πο 92 Δι γι 


$ 
1 


' 2 ' ' ee ' 2 
(A544? (Ai 1) (5i) 542 (Ai 44? 
ccu a ρω eee 
1151 111 itl 


1:63 


2 2 











πα “i Ai og M 
64-5; 1 IH SOM DOE Tal 
= 2 
μι ο -υ 
q E (I-12) 
111 
b. At level i = 1: 
λα κ [UMP 230000 
΄ 192 96 192 
= x _ 2 στι 
io). «Se h a 
i 64 192 64-Δ. 
c At level i - 2 
Uu, «A! ον. ΠΠ. (A*) CA?) 
T(2,1) = 1 2, a ρε, E 2 is 3.1» 
j 64 2 192 192°A, 24.54 
u,- (A5)? 
ιαπ. (1-14) 
ΠΝ Tae ean ΤΙ 
m2 2) = l 2 Pu M emen qo om - m 3 2 
í 192 ΠΡ 192 4874. 
πο. Ie T 
too x — T sock 
3 3 3 
E Geet 
D αἱ 5 
T dE s (1-15. 
u, (a αὐ’ αὐ. ΝΡ 
να. (64-84 ` + ο a 12.8. 3 64-4. (1. 
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d. At level i = N-1: 


E 2 2 

u (A! .) - BS T) 

T(N+1,N) = ao --- N+1 
N+1 
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r E EL 
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64 192 
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APPENDIX J 


FORCED VERTICAL VECOC IHLY 


l. The contribution to the surface geopotential from the 
forced vertical velocity is GSi 


eel E Saa 
Ou sin (ππτ) Sin (ux) tices 
OM Sin (ux) t E 
° ! ; 22 
where ou is mountain geopotential (m/s ), t is time, and 


T 1s the total time to burlid the mountain. 1s a constant, 


ÓM 


που (3-2) 


where g is gravity and H, is the height of the mountain. 


M 
Hu 1s a parameter specified in each model. Hy = 750 meters 
in the thesis experiments. 
2. The time rate of change of oe 1S separated into sine and 


cosine components for use in the surface geopotential forecast 


equations, 


4—- = MTS, (t)°cos(ux) + MTS, (t) -sin (px) : (J-3) 
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Cio -— ου. X š | 
where Hm or + U.tc' 3x Equation (J-l) is substituted into 


equation (J-3) and the resultant expression is separated 


into sine and cosine equations. The equations to calculate 


the terms MTS, and MTS, are 


= PPP ere 
u z Og u sin (s=) CNET 

MTS, (t) = 7 (J-4) 
Usfc "MU E un 

and 
το 
Μ πο TE 

— sin (s=) cos (s>) COSE 

MTS, (t) = (J-5) 
0 Ἐπ 


These terms are calculated for each time step in the model's 


forecast subroutine. 
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APPENDIX K 


DIABATIC HEATING TERMS 


For the diabatic heating function defined in equation 
(3-1), the maximum heating occurs at Z = 0.50, the midpoint 
of the heated layer. Temperature and diabatic heating are 
defined for grid A at the staggered levels, while they are 
defined at the unstaggered levels for grids B and C. Conse- 
quently, the rate of heating differs between the staggered 
and unstaggered models. In these experiments, the heated 
layer is between Z = 0.40 and Z = 0.60, the heating rate is 
5.0 K/day, and only the cosine term, Ql, 1S nonzero in the 
heated layer. The value of the heating term is listed below 
for staggered and unstaggered levels for six, twelve and 


sixty-layer models. 


Grid A 60-Layer Models Grid B and C 60-Layer Models 
z Q1 z Q1 
02592 0.0 0.400 0.0 
0.408 0.986E-06 0.417 0 SO OUS 
0.425 0.847E-05 0.433 0.145E-04 
0.442 0.214E-04 0.450 0.289E-04 
0.458 0.364E-04 0.467 0.434E-04 
0.475 0.419 3E-04 0.483 0.540E-04 
0.492 0.569E-04 0.500 0.579E-04 
0.508 0.569E-04 0.517 0.540E-04 
05525 0.493E-04 0:533 0.434E-04 
07522 0.364E-04 055550 0.289E-04 
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0.558 
0.575 
6: 592 
0.608 


1105512 L3aver Models 


E 

DEUS 
0.458 
0.542 


0.625 


Grid A 6-Layer Models 


p. 
92:50 
0.417 
0.583 
0750 


0.214E-04 

0.847E-05 

0.986E-06 
OSO 


Ql 

0.0 
0.364E-04 
0.364E-04 

0.0 


Ql 

0.0 
0.388E-05 
0.388E-05 

0.0 


0.567 0 1455=-04 
0.588 ο ος ο” 
0.600 Ὁ 


Grid B and C 12-Layer Models 


z. Q1 
0 535 0.0 
0.417 0.388E-05 
0.500 0.579 E-04 
0.585 0.388E-05 
0.667 0.0 


Grid B and C 6-Layer Models 


2 Ql 
0.333 0.0 
0.500 0.579E-04 
0.667 0.0 
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